subroutine def_p_conv(Pobs,x,xCell,s,n,p,rho,u,c,Pinf,dpdt,drdt,dudt,ttime,dtout,dt,m)
use Data_Cfoam, only: UFON
    implicit none

    real *8 Pobs(1),x(3),xcell(3),s,n(3),c,Pinf,Rhoinf,P,Rho,u(3),dPdt, drdt,dudt(3),ttime,dtout,dt,xdiff,ydiff,zdiff,rdiff
    real *8 Tobs0
    REAL *8 pi,Tobs,p01,p02,p03,p04,p05,p06,p07,p08,VflowMod,Mach,ER(3),ERmod,VflowR,VflowM,MachR,Mfactor,Lr,Lm,Vflow(3),Rtot(3)
    REAL *8 CosVflowR,SinVflowR
    REAL *8 RhoUn,RhoUnDot,LrDot,factor
!    real *8  UFON1 ! SK: + ONE VARIABLE FOR THE BACKGROUND FLOW

    integer,intent(out) :: m
    Vflow=0d0
    Vflow(1)=UFON
    Rhoinf=1.4*Pinf/c**2
    pi=3.141592653589793
    
    
    xdiff=x(1)-xCell(1)
    ydiff=x(2)-xCell(2)
    zdiff=x(3)-xCell(3)
    rdiff=sqrt(xdiff**2+ydiff**2+zdiff**2)
      
        
            
              VflowMod=(Vflow(1)**2+Vflow(2)**2+Vflow(3)**2)**0.5
              Mach=VflowMod/c
              
              if (VflowMod.eq.0) then
              CosVflowR=xdiff/rdiff
              SinVflowR=dsqrt(1-CosVflowR**2)
              else
              CosVflowR=(Vflow(1)*xdiff+Vflow(2)*ydiff+Vflow(3)*zdiff)/(rdiff*VflowMod)
              SinVflowR=dsqrt(1-CosVflowR**2)
              end if 
              factor=1d0/(1d0-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
              Tobs=ttime+rdiff/c*factor
    m=int(Tobs/dtOut)+1
! 1 -[0,dtout]
! 2 -[dtout,2*dtout]
     Tobs0=(m-1)*dtout
                           
!              m=int(Tobs/dtout) ! Tobs which we will use is Tobs(m)=m*dtout<=t3+rdiff/c
                             
              ER(1)=xdiff-Vflow(1)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
              ER(2)=ydiff-Vflow(2)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
              ER(3)=zdiff-Vflow(3)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
              
              ERmod=dsqrt(ER(1)**2+ER(2)**2+ER(3)**2)
              
              VflowR=(ER(1)*Vflow(1)+ER(2)*Vflow(2)+ER(3)*Vflow(3))/ERmod
              VflowM=VflowMod**2/c
              
              MachR=VflowR/c
              Mfactor=1+MachR

              !Lr=1/ERmod*((P-Pinf)*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3)) &
              !+Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3)))
              !Lm=1/c*((P-Pinf)*(Vflow(1)*n(1)+Vflow(2)*n(2)+Vflow(3)*n(3)) &
              !+Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(Vflow(1)*u(1)+Vflow(2)*u(2)+Vflow(3)*u(3)))
              
	      Lr=1/ERmod*((P-Pinf)*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3))&
	      +Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*(u(1)-Vflow(1))+ER(2)*(u(2)-Vflow(2))+ER(3)*(u(3)-Vflow(3))))
              Lm=1/c*((P-Pinf)*(Vflow(1)*n(1)+Vflow(2)*n(2)+Vflow(3)*n(3))&
	      +Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*&
	      (Vflow(1)*(u(1)-Vflow(1))+Vflow(2)*(u(2)-Vflow(2))+Vflow(3)*(u(3)-Vflow(3))))
      
             
              RhoUn=Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))
              RhoUnDot=drdt*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3)) &
              +Rho*(dudt(1)*n(1)+dudt(2)*n(2)+dudt(3)*n(3))
              !LrDot=1/ERmod*(dPdt*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3)) & 
              !+drdt*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3)) &
              !+Rho*(dudt(1)*n(1)+dudt(2)*n(2)+dudt(3)*n(3)) &
              !*(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3))+ &
	      !Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*dudt(1)+ER(2)*dudt(2)+ER(3)*dudt(3)))

	      LrDot=1/ERmod*(dPdt*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3))&
	      +drdt*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))&
	      *(ER(1)*(u(1)-Vflow(1))+ER(2)*(u(2)-Vflow(2))+ER(3)*(u(3)-Vflow(3)))&
	      +Rho*(dudt(1)*n(1)+dudt(2)*n(2)+dudt(3)*n(3))&
	      *(ER(1)*(u(1)-Vflow(1))+ER(2)*(u(2)-Vflow(2))+ER(3)*(u(3)-Vflow(3)))&
	      +Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*dudt(1)+ER(2)*dudt(2)+ER(3)*dudt(3)))
              
              p01=s*(Lr+Lm)/(ERmod**2*Mfactor**2)
              p02=s*LrDot/(c*ERmod*Mfactor**2)
              p03=-s*(Lr*c*MachR*Mfactor)/(c*ERmod**2*Mfactor**3)
              p04=-s*(RhoUn*(VflowR-VflowM))/(ERmod**2*Mfactor**2)
              p05=-s*(VflowR*RhoUnDot)/(c*ERmod*Mfactor**2)
              p06=s*(VflowR*RhoUn*c*MachR*Mfactor)/(c*ERmod**2*Mfactor**3)
              p07=s*(RhoUnDot)/(ERmod*Mfactor)
              p08=-s*(RhoUn*c*MachR*Mfactor)/(ERmod**2*Mfactor**2)
            
     if (Tobs-dt.lt.Tobs0) then
!  Pobs(m-1)=Pobs(m-1)+(p1+p2+p3+p4+p5)*(Tobs0-(Tobs-dt))/dtout
!   Pobs(m)=Pobs(m)+(p1+p2+p3+p4+p5)*(Tobs-Tobs0)/dtout
  Pobs(m-1)=Pobs(m-1)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*(Tobs0-(Tobs-dt))/dtout
  Pobs(m)=Pobs(m)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*(Tobs-Tobs0)/dtout
     else
              Pobs(m)=Pobs(m)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*dt/dtout
!    Pobs(m)=Pobs(m)+(p1+p2+p3+p4+p5)*dt/dtout
     end if
        
    return
end subroutine def_p_conv



!subroutine def_p_conv(Pobs,x,xCell,s,n,p,rho,u,c,Pinf,dpdt,drdt,dudt,ttime,dtout,dt,m)
!    implicit none
!
!    real *8 Pobs(1),x(3),xcell(3),s,n(3),c,Pinf,Rhoinf,P,Rho,u(3),dPdt, drdt,dudt(3),ttime,dtout,dt,xdiff,ydiff,zdiff,rdiff
!    real *8 Tobs0
!    REAL *8 pi,Tobs,p01,p02,p03,p04,p05,p06,p07,p08,VflowMod,Mach,ER(3),ERmod,VflowR,VflowM,MachR,Mfactor,Lr,Lm,Vflow(3),Rtot(3)
!    REAL *8 CosVflowR,SinVflowR
!    REAL *8 RhoUn,RhoUnDot,LrDot,factor
!    integer,intent(out) :: m
!    Vflow=0d0
!    Vflow(1)=0d0
!    Rhoinf=1.4*Pinf/c**2
!    pi=3.141592653589793
!    
!    
!    xdiff=x(1)-xCell(1)
!    ydiff=x(2)-xCell(2)
!    zdiff=x(3)-xCell(3)
!    rdiff=sqrt(xdiff**2+ydiff**2+zdiff**2)
!      
!        
!            
!              VflowMod=(Vflow(1)**2+Vflow(2)**2+Vflow(3)**2)**0.5
!              Mach=VflowMod/c
!              
!              if (VflowMod.eq.0) then
!              CosVflowR=1
!              SinVflowR=0
!              else
!              CosVflowR=(Vflow(1)*xdiff+Vflow(2)*ydiff+Vflow(3)*zdiff)/(rdiff*VflowMod)
!              SinVflowR=dsqrt(1-CosVflowR**2)
!              end if 
!              factor=1d0/(1d0-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
!              Tobs=ttime+rdiff/c*factor
!    m=int(Tobs/dtOut)+1
!! 1 -[0,dtout]
!! 2 -[dtout,2*dtout]
!     Tobs0=(m-1)*dtout
!                           
!!              m=int(Tobs/dtout) ! Tobs which we will use is Tobs(m)=m*dtout<=t3+rdiff/c
!                             
!              ER(1)=xdiff-Vflow(1)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
!              ER(2)=ydiff-Vflow(2)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
!              ER(3)=zdiff-Vflow(3)*rdiff/c/(1-Mach**2)*((1-(Mach*SinVflowR)**2)**0.5-Mach*CosVflowR)
!              
!              ERmod=dsqrt(ER(1)**2+ER(2)**2+ER(3)**2)
!              
!              VflowR=(ER(1)*Vflow(1)+ER(2)*Vflow(2)+ER(3)*Vflow(3))/ERmod
!              VflowM=VflowMod**2/c
!              
!              MachR=VflowR/c
!              Mfactor=1+MachR
!
!              Lr=1/ERmod*((P-Pinf)*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3)) &
!              +Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3)))
!              Lm=1/c*((P-Pinf)*(Vflow(1)*n(1)+Vflow(2)*n(2)+Vflow(3)*n(3)) &
!              +Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(Vflow(1)*u(1)+Vflow(2)*u(2)+Vflow(3)*u(3)))
!                            
!              RhoUn=Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))
!              RhoUnDot=drdt*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3)) &
!              +Rho*(dudt(1)*n(1)+dudt(2)*n(2)+dudt(3)*n(3))
!              LrDot=1/ERmod*(dPdt*(ER(1)*n(1)+ER(2)*n(2)+ER(3)*n(3)) & 
!              +drdt*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3)) &
!              +Rho*(dudt(1)*n(1)+dudt(2)*n(2)+dudt(3)*n(3)) &
!              *(ER(1)*u(1)+ER(2)*u(2)+ER(3)*u(3))+ &
!        Rho*(u(1)*n(1)+u(2)*n(2)+u(3)*n(3))*(ER(1)*dudt(1)+ER(2)*dudt(2)+ER(3)*dudt(3)))
!              p01=s*(Lr+Lm)/(ERmod**2*Mfactor**2)
!              p02=s*LrDot/(c*ERmod*Mfactor**2)
!              p03=-s*(Lr*c*MachR*Mfactor)/(c*ERmod**2*Mfactor**3)
!              p04=-s*(RhoUn*(VflowR-VflowM))/(ERmod**2*Mfactor**2)
!              p05=-s*(VflowR*RhoUnDot)/(c*ERmod*Mfactor**2)
!              p06=s*(VflowR*RhoUn*c*MachR*Mfactor)/(c*ERmod**2*Mfactor**3)
!              p07=s*(RhoUnDot)/(ERmod*Mfactor)
!              p08=-s*(RhoUn*c*MachR*Mfactor)/(ERmod**2*Mfactor**2)
!            
!     if (Tobs-dt.lt.Tobs0) then
!!  Pobs(m-1)=Pobs(m-1)+(p1+p2+p3+p4+p5)*(Tobs0-(Tobs-dt))/dtout
!!   Pobs(m)=Pobs(m)+(p1+p2+p3+p4+p5)*(Tobs-Tobs0)/dtout
!  Pobs(m-1)=Pobs(m-1)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*(Tobs0-(Tobs-dt))/dtout
!  Pobs(m)=Pobs(m)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*(Tobs-Tobs0)/dtout
!     else
!              Pobs(m)=Pobs(m)+1/(4*pi)*(p01+p02+p03+p04+p05+p06+p07+p08)*dt/dtout
!!    Pobs(m)=Pobs(m)+(p1+p2+p3+p4+p5)*dt/dtout
!     end if
!        
!    return
!end subroutine def_p_conv
